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The basic formulation describing quadratic mode coupling in rotating Newtonian stars is pre¬ 
sented, focusing on polar modes. Due to the Chandrasekhar-Friedman-Schutz mechanism, the 
/-mode (fundamental oscillation) is driven unstable by the emission of gravitational waves. If the 
star falls inside the so-called instability window, the mode’s amplitude grows exponentially, until it 
is halted by nonlinear effects. Quadratic perturbations form three-mode networks inside the star, 
which evolve as coupled oscillators, exchanging energy. Coupling of the unstable /-mode to other 
(stable) modes can lead to a parametric resonance and the subsequent saturation of its amplitude, 
thus suppressing the instability. The saturation point determines the amplitude of the gravitational- 
wave signal obtained from an individual source, as well as the evolutionary path of the latter inside 
the instability window. 
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I. INTRODUCTION 

With a mass of the order of the solar mass and a radius 
of about 10 km, neutron stars constitute nature’s high- 
energy laboratories, from which the behavior of matter at 
such extreme conditions could be deduced. The neutron 
star equation of state is yet to be determined and remains 
one of the most significant questions in astrophysics. The 
serendipitous discovery of the first pulsar by Hewish and 
Bell in 1967 signified the onset of neutron star astronomy, 
which has provided some constraints for the masses, the 
radii, and the rotation periods of neutron stars. 

Nevertheless, these observations are not enough to in¬ 
fer the equation of state. A method that could be used to 
further probe neutron stars is asteroseismology, namely, 
the study of stellar oscillations [1, 2]. Especially after 
the realization that stellar oscillations can be driven un¬ 
stable by the emission of gravitational radiation [3, 4], 
the field of gravitational wave asteroseismology was de¬ 
veloped rapidly; detection of gravitational waves from 
nonradial stellar oscillations could provide information 
about the neutron star interior [5-9]. 

The Chandrasekhar-Friedman-Schutz (CFS) instabil¬ 
ity, however, grows on long time scales and, to make 
things worse, it is suppressed by viscosity [10, 11]. For 
the /-modes, which are the fundamental oscillations of 
the star and the best gravitational wave emitters, this 
leaves only a small portion of the parameter space, where 
the instability is active. In the late 1990s, it was realized 
that another class of oscillations, the r-modes, is unstable 
for a much larger parameter range [12-16]. The r-modes 
are related to horizontal motions of the fluid, much like 
Rossby waves in the Earth’s atmosphere and oceans, and 
exist only in rotating stars [17]. Moreover, they have 
shorter growth times, compared to the /-modes. As a 
result, the r-mode instability was considered as the most 
promising gravitational wave source. 

Consequent studies on the r-mode instability naturally 
raised the question of the maximum amplitude that the 


oscillation can attain, before it is halted by nonlinear 
effects. Coupling of the unstable r-mode to other modes 
of the star can work as an energy drain and saturate 
the instability. The results of these studies were quite 
disappointing, from a gravitational-wave-detection point 
of view: the r-mode saturation amplitude is, in fact, quite 
lower than expected, or, at least, hoped [18-21]. 

Determining the saturation amplitude of the unsta¬ 
ble r-mode is also important for neutron star evolution. 
Whether the star is newborn or a member of a low-mass 
x-ray binary system (LMXB), its evolution depends on 
the value of the saturation amplitude [22, 23]. When 
the star enters the instability region, it loses angular mo¬ 
mentum, due to gravitational wave emission, which could 
possibly explain the upper limit in the observed neu¬ 
tron star rotational frequencies [14, 24-27] (about 700 
Hz [28]). 

Even though the r-mode instability is active in a much 
larger part of the parameter space, the /-mode instability 
could still be significant, especially for newborn neutron 
stars. Furthermore, the fact that the r-mode saturation 
amplitude is not expected to be high renders the study of 
the /-mode quite important: if the /-mode is not satu¬ 
rated at such low amplitudes, then it could be a possible 
gravitational wave source and, thus, provide much infor¬ 
mation about the neutron star equation of state. Up until 
now, the evolution of the /-mode instability in the non¬ 
linear regime has been performed only via hydrodynamic 
simulations [29-31]. However, since the growth time of 
the instability is, in general, quite long, it is very hard 
for nonlinear simulations to track the mode evolution for 
such a long time. 

Recent work [32] suggests that, should the /-mode sat¬ 
urate at reasonably high amplitudes, the gravitational 
wave signal from a source in the Virgo cluster, under¬ 
going the /-mode instability, could be detectable by the 
Einstein Telescope. A more promising source is related 
to supramassive configurations (exceeding the maximum 
mass of a nonrotating star), which could be the outcome 


2 


of a neutron star merger. Such stars would be stable 
only for rotation rates close to the Kepler limit (mass- 
shedding limit). The /-mode instability is expected to 
grow really quickly in these objects and the gravitational 
wave signal could even reach the sensitivity of Advanced 
LIGO, with a quite promising event rate [33]. 


As opposed to the r-mode, where the oscillation com¬ 
prises horizontal fluid motions, the /-mode is dominated 
by a radial component and large-scale density varia¬ 
tions, which makes it a more efficient gravitational wave 
emitter. However, the so-called instability window is 
much smaller for the /-mode. This is a region in the 
“temperature-rotation rate” plane, where the instabil¬ 
ity is not suppressed by viscous effects. By expanding 
a perturbation in its multipole moments (described by 
the spherical harmonics T]™) we see that higher multi¬ 
poles become unstable at lower rotation rates. On the 
other hand, lower multipoles emit gravitational waves 
more efficiently, but might not become unstable at all. 
The instability window of the I = m = 2 /-mode is sig¬ 
nificant only for models with quite stiff equations of state, 
whereas I = m = 3 and 4 /-modes have larger windows, 
but might not grow very fast. 


Applying the same methodology as for the r-mode, we 
can determine the amplitude at which the /-mode insta¬ 
bility is saturated by nonlinear effects. This work has 
been divided into two parts. In the first part, included in 
the present paper, we will present the theoretical frame¬ 
work of the problem. Its application to various stellar 
models will be presented in a subsequent paper. 


The paper is organized as follows: in Sec. II we present 
the formalism that gives rise to the various oscillation 
modes in the star, using linear perturbations. We dis¬ 
cuss the method with which one can acquire the oscilla¬ 
tion spectrum in the nonrotating limit, and then we add 
rotation in a perturbative way (slow-rotation approxima¬ 
tion) and present its main implications. In Sec. Ill we 
give a short overview of the CFS instability and how it is 
manifested in the /-mode. In Sec. IV we review the for¬ 
malism which describes quadratic perturbations and de¬ 
rive the conditions under which coupled-mode networks 
can arise. Furthermore, these networks are subjected to 
a stability analysis, which determines whether saturation 
can be achieved by the system or not. Derivations of sev¬ 
eral formulas in this section are addressed in Appendices. 
In Appendix A we derive the equations of motion, includ¬ 
ing quadratic perturbations, whereas in Appendix B we 
give the expression for the three-mode coupling coeffi¬ 
cient. Appendix C contains a study of a coupled three¬ 
mode network, using the multiscale method, as well as 
the details of the stability analysis mentioned above. Fi¬ 
nally, Sec. V concludes the paper with some discussion. 


II. THE OSCILLATION MODES—LINEAR 
PERTURBATION SCHEME 

Stellar oscillation modes can be divided in two gen¬ 
eral categories: polar (or spheroidal) modes and axial 
(or toroidal) modes. Expanding the displacement vec¬ 
tor field of an arbitrary perturbation in vector spherical 
harmonics, we get 

i 

i{r,e,(k) = Y. E Wr{r)Yr{0,^)er ( 2 . 1 ) 

I m——l 

+ Vr{r)VYr{e, /.) +ur{r)er x </')], 

where (r, 0 , />) are the spherical polar coordinates, 
(ef,,eg,e 0 ) is the orthonormal basis, and 1 )™ are the 
spherical harmonics. Then 

• polar modes: UJ^ = 0 

' as D ^ 0 , 

• axial modes: Vj™ = IF/” = 0 

D being the stellar rotation rate, /-modes, as well as p- 
(acoustic waves) and g-modes (gravity waves), are exam¬ 
ples of polar modes. They constitute the “regular” mode 
spectrum of a star and have finite frequencies in the non¬ 
rotating limit, r-modes, on the other hand, are axial and 
become trivial in the nonrotating limit, where their fre¬ 
quencies vanish (for a detailed presentation of oscillation 
modes, cf. for instance, Refs. [1, 2]). The picture above 
slightly changes in the case of zero-buoyancy stars, g- 
modes, which are caused by the presence of buoyancy, 
become trivial too. The result of this “mixture” of triv¬ 
ial modes (r- and 5 -modes) is another class of modes, 
called hybrid modes, which have both polar and axial 
components in the nonrotating limit. In the special case 
where I = m one obtains the “classical” r-modes, which 
are purely axial [34]. 

Assuming a star which is uniformly rotating with an 
angular velocity D, the fluid equations, in the frame ro¬ 
tating with the star, are 

+ V • (pv) = 0, (2.2) 

dv Vd 

— —h ("W • V)-?; 2fl X V fl X (fl X r) = - 

at p 

(2.3) 

and 

= 47rGp, (2.4) 

where p is the density, p the pressure, v the velocity, 
the gravitational potential and G the gravitational con¬ 
stant. The system above has to be supplemented with 
an equation of state p = p{p,p), where p usually cor¬ 
responds to entropy or composition and depends on the 
density. By considering “small” perturbations imposed 
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on the equilibrium state, these equations are written as 

(2.5) 


^ + V • (pSv) = 0, 


dSv VSp 

— -h 2fl X Sv = - 

ot p 


^Sp-VS^, 


= AirGSp, 


( 2 . 6 ) 

(2.7) 


where and are the radial and horizontal components 
of respectively. It should be noted that, since operator 
C is Hermitian [35], the solutions to Eq. (2.13) (with 
vanishing B) are orthogonal, i.e. 

= J pC- (2.15) 


and 


where 


Ap ^ Ap / i91np \ A/i 

p ^ p \d\ap) p p ’ 


El = 


/ 91np\ 
\d\np) 


( 2 . 8 ) 


(2.9) 


In the equations above 5 denotes a Eulerian perturbation 
and A corresponds to a Lagrangian perturbation. The 
former monitors changes in a particular point in space, 
whereas the latter refers to changes in a given fluid ele¬ 
ment. The two are related by A/ = 5f + • V)/, where 

^ is the Lagrangian displacement of the fluid element 
[1, 35]. 

By definition, Av = d^/dt = ^ -I- (w • V)^, but, since 

= 0 in the background, At> = ^ = 5v. Then, the 
perturbed Euler equation (2.6) can be written as [35] 



| + B(|)+C(0 = 0, 

(2.10) 

where 


B{i) = 2nxt 

(2.11) 

and 



(2.12) 


Operator C can be written in terms of ^ by using 
Eqs. (2.5), (2.7), and (2.8) to replace the perturbations 
dp, dd), and dp, respectively (cf. for example. Sec. II B 
in Ref. [18], or Sec. 2.1 in Ref. [35]). 

Seeking solutions of the form ^{r,t) = ^(r)e“*, where 
ui denotes the frequency of a mode in the corotating 
frame, Eq. (2.10) is written as 

+ = 0. (2.13) 


This is the eigenvalue equation which needs to be solved, 
supplemented with the appropriate boundary conditions, 
in order to obtain the mode spectrum of the star. 


A. The nonrotating limit 

Equation (2.13) is simplified significantly in the ab¬ 
sence of rotation, since operator B vanishes. Then, ac¬ 
cording to Eq. (2.1), the displacement vector of a polar 
mode is 

(2.14) 


where the indices in ^ denote different solutions, d^p is 
the Kronecker delta, and the star denotes complex con¬ 
jugation. Since all perturbative quantities are functions 
of I, they can all be expressed as 

d/(r,0,</>,t)=d/(r)yr(^,</')e“‘- 

Hence, a separation of variables is possible and the prob¬ 
lem is reduced to calculating the radial dependence of 
the perturbation [ 1 ]. 

A sample from the polar mode spectrum of a poly¬ 
tropic star is presented in Fig. 1. Each mode is generally 
described by three numbers: its overtone n, its degree 
I, and its order m. When rotation is absent, the mode 
frequencies do not depend on m (see Sec. IIB). The /- 
mode (n = 0 ) lies between its overtones (n > 0 ), the 
high-frequency p-modes and the low-frequency p-modes. 
p-modes are pushed towards zero as the effects of buoy¬ 
ancy become less and less important, until they finally 
vanish for zero-buoyancy stars. Departure from the zero- 
buoyancy case can be a result of stratification (composi¬ 
tion gradients) or deviations from isentropy (star with a 
finite temperature) [36]. 

This behavior can be described by the so-called 
Schwarzschild discriminant, which is given by 

^ dlnp 1 dlnp 

dr Ti dr ’ 

where Ti is the adiabatic exponent, defined in Eq. (2.9). 
If the star obeys a simple polytropic equation of state p = 



FIG. 1. Polar mode spectrum of a star obeying a poly¬ 
tropic equation of state with F = 2. The adiabatic expo¬ 
nent Fi is equal to 2.1. Mode frequencies, which scale as 
uj = ul-sjGMIB?, are plotted against the mode degree 1. 
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Kp^ (where K and T are constants), the Schwarzschild 
discriminant becomes 

^ Fi — rdlnp 

Then, if F = Fi the star exhibits no convective phenom¬ 
ena (zero-buoyancy case). On the other hand, F < Fi 
(F > Fi) denotes convective stability (instability), i.e. 
oscillatory (unstable) g-modes. 

If the equation of state is described by the more general 
relation p = p{p, p ), the occurrence of convective phe¬ 
nomena is parametrized through p. The condition for the 
existence of (;-modes is /S.p = 0 [cf. Eq. (2.8)]. If p cor¬ 
responds to the composition, this condition means that 
the composition of a displaced fluid element is “frozen”; 
weak interaction processes need more time than an os¬ 
cillation period to restore /3-equilibrium between the dis¬ 
placed fluid element and the surrounding matter. On 
the other hand, if p corresponds to entropy, it means 
that the fluid displacement occurs adiabatically. The 
Schwarzschild discriminant, as a function of p, is given 

by 

^ I /91np\ dlnp 
Fi V51n/rj^ dr 

B. The slow-rotation approximation 

Taking rotation into account, the situation changes sig¬ 
nificantly. The equilibrium configuration no longer ex¬ 
hibits spherical symmetry and an oscillation mode cannot 
be described by a single spherical harmonic. Typically, 
Eq. (2.13) has to be solved from scratch. However, ro¬ 
tation can also be introduced perturbatively, namely, by 
considering the effects of rotation to the various quan¬ 
tities as perturbations. Rotation affects polar modes in 
two ways. First, it lifts the (2Z-|-l)-fold degeneracy in the 
eigenfrequency of each mode, by introducing a Zeeman- 
like splitting. The eigenfrequency now depends on both 
the degree I and the order m, as opposed to the nonrotat¬ 
ing limit, where it is degenerate in m. Second, rotation 
distorts the equilibrium structure of the star, which also 
changes the mode frequencies. An additional effect of ro¬ 
tation is, as discussed before, the appearance of a whole 
different class of modes, the inertial modes (like the r- 
mode), whose restoring force is the Coriolis force. 

Mode splitting is already introduced as a first-order 
effect, whereas equilibrium distortion is a second-order 
effect. Higher-order effects also become important for 
large rotational velocities, but the analysis is quite cum¬ 
bersome even at second order in H. A third-order per¬ 
turbation formalism was developed in Ref. [37] , where an 
interesting case of near degeneracy was observed. Nev¬ 
ertheless, we stopped at quadratic perturbations in O, 
keeping in mind that higher-order effects could be signif¬ 
icant at the near-Kepler angular velocities that we are 
interested in. 


Eigenfrequencies, eigenfunctions, as well as equilib¬ 
rium quantities, are expanded as 

oj = ujq u;i(H) -1- cu 2(H^) -t- (Il(H^), 

^ = ^o + ^im + U^^) + o{n^), 
p = Po + P2{n^) + oin‘^). 

Substituting these in Eq. (2.13) and distinguishing be¬ 
tween first- and second-order terms, we obtain [18] 

~ -I- Co(^i) — 2ujouii^o + iujoBi{^o) = 0 (2.16) 

and 

“^0^2 + ^ 0 (^ 2 ) — 2a;oWi^i -I- — 2ujoU}2^o 

— ^1^0 + *wiBi(4o) + ^ 2 (^ 0 ) = 0, (2-17) 

respectively. From the above, we find the 0{ft) and 
O(H^) corrections to the eigenfrequencies. The first is 
rather simple and is given by 

wi = mCiH, (2.18) 

where 

^ ^ I[‘2^r^h + il]pr^dr 

^ I[& + Ki + ^)^l]pr^dr' 

The second is more complicated and has the general form 
[38] 

uj2=C2n^ = {X+ m^Y)n^, (2.19) 

where X and Y include corrections due to the distortion 
of the equilibrium and due to the effects of the Coriolis 
force. The effect of rotation on the mode eigenfrequencies 
(up to second order) can be seen in Fig. 2. 



FIG. 2. Eigenfrequency of the 1 = 2 /-mode (in the corotat¬ 
ing frame), as a function of the rotation rate H. Each line 
corresponds to a different value of m. As in Fig. 1, a poly¬ 
trope with F = 2 and Fi =2.1 was used. The mode frequency 
scales as cj = ujj y/GM/, whereas the rotational velocity is 
normalized to the Kepler limit Hk- 
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As for the eigenfunctions, rotation couples polar modes 
to axial modes, as well as other polar modes. This means 
that a mode cannot be any more described by a single 
spherical harmonic, which makes the situation more com¬ 
plicated. Since operator B is nonvanishing in this case, 
the solutions to Eq. (2.13) do not obey the orthogonality 
relation (2.15). Instead, they satisfy a modified orthogo¬ 
nality condition, given by^ [18] 

(oJa + W/3)(|a, - (^a, i®(l/3)) = baSafi. (2.20) 


III. THE /-MODE INSTABILITY 


As it was discovered by Chandrasekhar [3] and rig¬ 
orously proven by Friedman and Schutz [4], oscillation 
modes can be driven unstable by the emission of gravi¬ 
tational radiation, if the star is rotating rapidly enough. 
Every mode can be thought of as having a prograde (de¬ 
noted by — |to|) and retrograde (denoted by jmj) compo¬ 
nent. Should the star rotate sufficiently fast, it can drag 
the retrograde component towards the direction of rota¬ 
tion, making it appear as prograde to a distant observer. 
Emission of gravitational waves by the perturbation can 
then act as a driving mechanism, increasing the mode 
energy. This can be seen by the standard multipole ex¬ 
pansion of the power radiated in the form of gravitational 
waves (GW) [39] 


- {\SDrf + \6Jrf) . 

/ GW 

(3.1) 

As one can see, the power emitted is negative (gravita¬ 
tional radiation damps the mode), unless uj{ui — mQ) < 0, 
in which case the energy of the mode is increased. The 
onset of the instability occurs when uj/m = fl, namely 
when the pattern speed of the mode matches the angu¬ 
lar velocity of the star. The angular velocity at which 
this happens is usually called critical. Alternatively, 
w — toH can be thought of as the mode frequency, mea¬ 
sured in an inertial frame {uj is the corotating frame fre¬ 
quency). Then, the instability sets in at the point when 
the inertial-frame frequency changes sign. 

In Eq. (3.1), Ni is a constant given by 


47rG il + l){l + 2) 
l{l - 1) [{21 + l)Uf 


(3.2) 


(c being the speed of light), whereas and 5J™ de¬ 
note the mass and current multipole moments, respec¬ 
tively. The /-mode radiates mainly via the former,^ 


which are given by 

50-^ = j r^5pYi*'^d?r. (3.3) 

Finally, the lower limit of the sum is given by /min = 
max(2 , |to|). 

Depending on the equation of state, all the / = m 
/-modes can become unstable. However, various dissi¬ 
pation mechanisms are expected to act against the CFS 
instability. Responsible for the dissipation of the /-mode 
are mainly bulk and shear viscosity (BV and SV), and 
their contributions are given by [11] 

and 

respectively. Here, 6a'^^ is the stress tensor and is given, 
in terms of the velocity perturbations, by 

^ -f , (3.6) 

Sa = Vi6v\ (3.7) 

Pij being the spatial metric tensor. / and t] are the 
bulk and shear viscosity coefficients, which depend on the 
equation of state (cf. for instance. Ref. [40]). Bulk viscos¬ 
ity is a result of the fluid trying to restore /3-equilibrium 
and operates at high temperatures, as opposed to shear 
viscosity, which is due to particle scattering and is dom¬ 
inant at low temperatures. 

For normal nuclear matter, comprising (nonsuperfluid) 
neutrons, (nonsuperconducting) protons, and electrons, 
neutron collisions make the biggest contribution to shear 
viscosity, and the two coefficients are given by [11, 41, 42] 

( = 6 X gcm-^s"^ (3.8) 

and 

rj = M7gcm-^ s-\ (3.9) 

where T is the stellar temperature and all the quanti¬ 
ties have cgs units. For superfluid nuclear matter an¬ 
other dissipation mechanism dominates, called mutual 
friction. This is expected to occur for temperatures 
< lO^K and suppresses the instability very efficiently 
[43]. Here, we only consider normal nuclear matter; as 
shown by Ref. [32], the star may never enter the su¬ 
perfluid region, since neutrino cooling is balanced by 
the oscillation-induced viscous heating before the star 
reaches the transition temperature.^ 


^ Note that Ref. [18] uses a different ansatz for ^{r, t), i.e. t) = 
hence the sign difference in the second term. 

^ Current multipole moments become significant in the case of the 
r-modes (cf. for example, Ref. [14]). 


^ Reference [32] uses E = 10“^ MQP'B? for the saturation energy 
of the /-mode. However, if the saturation energy is smaller, 
viscous heating due to the oscillation balances neutrino cooling 
at lower temperatures. 
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FIG. 3. Instability windows of the I = m = 3 and I = m = i 
/-modes, for a polytropic model with F = 2 and Fi = 2.1 
(the I = m = 2 /-mode does not become unstable for this 
model). Fiducial values were used for the mass and radius of 
the star, i.e. M = 1.4 Mq and R = 10 km. The angular veloc¬ 
ity is normalized to the Kepler limit Hk. These curves were 
not produced using the slow-rotation formalism described in 
Sec. IIB, because the modes fail to become unstable in this 
approximation. Also, although this model does not favor the 
instability, making use of realistic equations of state and rel¬ 
ativity can push the windows to quite lower angular velocity 
values [44, 45]. 


The instability is active only if the total energy rate of 
the mode is positive, i.e. 


dE 

dt 


V dt / \ dt y gy 



(3.10) 


path of the neutron star inside the instability window. 
After the star enters the instability window, it cools 
down, until neutrino cooling is balanced by viscous heat¬ 
ing due to the oscillation. Then, it descends the instabil¬ 
ity window at almost constant temperature, by losing an¬ 
gular momentum. However, magnetic braking also slows 
down the star, competing with gravitational radiation; as 
shown in Ref. [32], the instability may not have enough 
time to grow, if the spin-down of the star is dominated 
by the magnetic torque. 

As in previous work for the r-mode instability [18- 
21, 23], we will consider quadratic perturbations and 
study their effects in the evolution of the /-mode. Even 
higher than second-order terms could, in principle, be im¬ 
portant at large oscillation amplitudes, but the complex¬ 
ity of the formulation and the requirements of our prob¬ 
lem allow us to choose simplicity over accuracy. Work 
that also includes cubic nonlinearities can be found in 
Refs. [46, 47]. Also, for a more general investigation 
of systems with quadratic and cubic nonlinearities the 
reader is referred to Chapter 6 of Ref. [48]. 


A. Mode decomposition 

As mentioned in Sec. II A, operator C of Eq. (2.13) 
is Hermitian. This means that, in the nonrotating limit 
(where B vanishes), any perturbation, described by the 
displacement vector can be decomposed as 

9a (4-1) 


By solving this inequality, one obtains the instability win¬ 
dow of the mode, namely the region in the T-il plane 
where the mode is CFS unstable (Eig. 3). Once the star 
enters this area, the amplitude of the mode will grow, un¬ 
til such a point where nonlinear effects become important 
and saturate it. This will be discussed in the following 
section. 


IV. MODE COUPLING—QUADRATIC 
PERTURBATION SCHEME 


where ^a(r) is a solution to Eq. (2.13) (with vanishing B) 
and represents the eigenfunction of an oscillation mode, 
whereas qait) is the amplitude coefficient. In the case of 
polar modes, this eigenfunction is given by Eq. (2.14). 

If rotation is included, operator B is nonvanishing and 
the solutions to Eq. (2.13) are not orthogonal, in general. 
However, instead of a configuration space mode expan¬ 
sion, like Eq. (4.1), one can use a phase space mode ex¬ 
pansion [49]. Then, a perturbation can be decomposed 
as [18] 


Considering the perturbations as small, the modes of 
the star are uncoupled oscillations (in the nonrotating 
limit). This is a result of the linear approximation used 
to define them (cf. Sec. H). However, as the amplitude of 
the unstable mode grows, the linear approximation fails 
to accurately describe it; higher-order terms are bound 
to play an important role in the amplitude evolution, 
since they introduce mode coupling. The result of this 
interaction of the unstable mode with other modes is the 
eventual saturation of the unstable mode’s amplitude. 

The actual value of this saturation amplitude is mainly 
important for two reasons. First, it sets the maximum 
amplitude of the gravitational wave signal obtained from 
the unstable mode. Second, it affects the evolutionary 




a ^ 

+ Q*Jt) 


ia{r) 

_ iWala(r) 

Cir) 





(4.2) 


This result was obtained by using the fact that both 
(oJcia) and {—oJaffia) are solutions to Eq. (2.13), as 
well as assuming that ${r,t) is real. 


B. Equations of motion 

Including second-order perturbative terms in 
Eq. (2.10), one obtains the quadratic equation of 
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motion, which can be generally written as 


amplitude equations of motion for the three modes are 


i' + B(O+C(O+Ar=0, (4.3) 

where A/” collectively denotes all 0{^^) terms. Substitut¬ 
ing Eq. (4.2), and using the eigenvalue equation (2.13) 
and the orthogonality condition ( 2 . 20 ), we get 


Qait) = f (4.4) 

This is the equation of motion for the amplitude of the 
mode Qa- If quadratic terms are ignored (or, equiva¬ 
lently, if the perturbation is small), then the amplitude 
Qa is constant, since there is no interaction with other 
modes. However, a nonzero Af couples the mode denoted 
by a with other modes, leading to an energy exchange be¬ 
tween them. For a derivation of the equations of motion 
(4.3) and (4.4), cf. Appendix A. 

By further replacing Eq. (4.2) in A/", we obtain 


Qa(i) 


fEE[ ‘E ctjS'yQ^ ^ 

“ 7 


(4.5) 


where J- denotes the coupling coefficient and is generally 
given by 


•^a/37 = (la,A/'(|/3,|T,)). (4.6) 

Borrowing the notation of Ref. [18], a bar over an index 
means that the corresponding mode eigenfunction in Af 
has to be complex conjugated and its frequency sign re¬ 
versed. The explicit form of the coupling coefficient is 
given in Appendix B. 

Observing Eq. (4.5), we see that modes couple in 
triplets, which is a natural consequence of the quadratic- 
perturbation approximation. This does not, however, re¬ 
strict the number of couplings for a single mode; if a 
mode couples to a pair of other modes, it can simultane¬ 
ously couple to other pairs as well. Also, one can notice 
that not all terms of Eq. (4.5) are equally significant. 
Rapidly varying terms do not contribute much on long¬ 
term dynamics and average to zero, as opposed to slowly 
oscillating components (this is proven by means of the 
multiscale method in Appendix Cl). Hence, couplings 
which really affect the mode amplitude evolution ought 
to satisfy a resonance condition, e.g. 

U>a = W/3 + W 7 + Aw, (4-7) 

where Aw is a small detuning (Aw Wi). Assuming 
such a relation between the mode frequencies, we can 
single out a mode triplet and follow its evolution. The 


Qa = (4.8a) 

Dq, 

Qp = (4.8b) 

bp 

Q 7 = (4.8c) 

So far, we have assumed that the modes are simply har¬ 
monic oscillations, unaffected by any growth/damping 
mechanisms. However, as discussed in the previous sec¬ 
tion, all the modes are influenced by various effects, such 
as gravitational radiation and viscosity. The majority of 
the modes is damped by these mechanisms, whereas a 
handful of modes can become unstable and grow, for a 
certain parameter range. 

Such effects are often parametrized by the imaginary 
part of the oscillation frequency. But we have hitherto 
assumed that mode frequencies are real, since no such 
effects were introduced in our equations. So, in order to 
calculate growth/damping rates, we will use the defini¬ 
tion of the corotating-frame mode energy, which is given 
by [18] 

— [gal ^aba 

= IQal^UJa [2w„(|„,la) " fB(^„))] . (4.9) 


This is a quadratic functional of 4, so, if 7 is the imagi¬ 
nary part of the frequency, then 

^ = 2jaEa. (4.10) 

at 

Formulas for dE /dt for the various mechanisms were pro¬ 
vided in the previous section, so we can calculate the 
growth/damping rate 7 for a particular mode. 

Incorporating the growth/damping rates in Eqs. (4.8), 
we get 

jOJ 

Qo. = iMo. + 

Dq; 

jqj 

Qp = ipQp + -^g;gc.e*^“‘, 

jqj 

where we also replaced the coupling coefficients with H = 
Eapj = Ep^a = E^ap (cf- Appendix B). 

Such three-mode systems can give an estimate of the 
effects of nonlinear coupling to the amplitude of an unsta¬ 
ble mode, like the /-mode. Such a mode, which we shall 
call “parent”, has 7 > 0 and has to be coupled to two 
“daughter” modes, which are linearly damped (7 < 0 ). 
The efficiency of the coupling depends on the value of 
the coupling coefficient "H, as well as on how close to 
resonance the three modes are. As we will see, some ad¬ 
ditional conditions have to be met, in order for the triplet 
to reach an equilibrium and saturate. 


(4.11a) 

(4.11b) 

(4.11c) 







(4.16) 


C. Mode normalization 


and 


For the amplitude coefficients of the modes Q to be 
meaningful, we first have to normalize all the modes ac¬ 
cording to some convention. By doing this, we will be 
able to compare the modes, using the same standards. 
The most popular normalization choice is to fix the mode 
energy (4.9) at unit amplitude to some arbitrary value 
Sunit, namely, 


for all modes. References [19-21] use i?unit = 
whereas Ref. [32] also uses Flunit = Mc^. The conver¬ 
sion between two different normalization choices can be 
straightforwardly written as 

IQal'i^unit = IQaPKnif (4-13) 

Using a normalization choice of the form (4.12), we can 
rewrite Eqs. (4.11) as 

Qa = laQa + (4.14a) 

Qp = JpQp + (4.14b) 

-^unit 

Q'Y = 77 Q 7 + . (4.14c) 

tiunit 

From this form of the amplitude equations of motion it 
is easier to see that the coupling coefficient % has units 
of energy. For the sake of generalization, though, we will 
be using Eqs. (4.11) in the subsequent sections.^ 


D. Coupling selection rules 

As we already mentioned, the three modes forming 
the coupled network have to obey a resonance condition, 
given by Eq. (4.7). The structure of the coupling coeffi¬ 
cient imposes two more conditions, which have to be met 
in order for the coupling to occur. 

As shown in Appendix B, the angular dependence of 
the zeroth-order component of the coupling coefficient 
has the form 

J j sin 6 >d 6 »d^, 

where Y)'" is the spherical-harmonic angular dependence 
of each mode [cf. Eq. (2.14)]. This integral is propor¬ 
tional to the Clebsch-Gordan coefficients (cf. for in¬ 
stance, Ref. [50]) and is nonzero if 

ma=mp + mj (4-15) 


— ^7 T 2A, 


where 


^ ^ 


and A — 0,1,... Aniax 



Equations (4.7), (4.15), and (4.16) constitute the selec¬ 
tion rules which the coupled mode triplet has to satisfy 
and restrict the search for possible couplings.® 


E. Parametric resonance instability 


As mentioned before, we are particularly interested in 
the case where an unstable parent mode ( 7 ^ > 0 ) is 
coupled to two damped daughter modes ( 7 / 3,7 < 0). In 
the beginning of the evolution, when the amplitudes are 
small, linear terms dominate: the amplitude of the par¬ 
ent grows and the amplitudes of the daughters decrease. 
At some point, nonlinear terms catch up and the parent 
starts pumping energy into the daughters. This point oc¬ 
curs when the parent exceeds a certain amplitude, called 
the parametric instability threshold. Such an interaction 
between the modes is an example of a parametric res¬ 
onance instability, i.e. an instability which can occur 
when the parameters of an oscillator vary in time (cf. for 
example. Ref. [51]). 

In order to obtain the parametric instability threshold, 
we take the daughters’ equations of motion (4.11b) and 
(4.11c) and ask what the value of the parent’s amplitude 
Qa should be, in order for the daughters’ amplitudes Q/ 3,7 
to start growing. Setting Qp,-y = Q/ 3,7 exp(jAa;t/2) and 
writing these equations in matrix form, we get [52] 


( *3/3 \ / 7/3 - i^oj/2 iQa'H/bp \ / Q/3 

V q; y V -iQiH/b^ 7 ^+ i^u/2 J q; 


{Qa is considered an unknown constant). The eigenval¬ 
ues of the system matrix are 


Ai ,2 - 2 


/ 9 47^2 

7/3-b 77 ± W (77 - 7/3-b iAw) +rirl* 3 c'^ 


6/3 b^ 


Then, for the system to admit a growing exponential so¬ 
lution, i.e. for the daughter modes to grow, the condition 
Re(A) > 0 has to be satisfied, for at least one of the eigen¬ 
values. This gives 


\Qo 


> 


7 / 3776/367 

772 


1 + 


Aw 


7/3+77 


(4.17) 


If one chooses a normalization of the form (4.12), they can simply 
replace with (^a'H/Eunit in t-h® following sections. 


^ It should be noted that, even though we evaluated the coupling 
coefficient in the nonrotating limit in Appendix B, these selection 
rules are valid to all orders in as shown by Ref. [18]. 
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which is the expression for the parametric instability 
threshold (PIT), i.e. the amplitude that the parent has 
to surpass so that the daughters will start growing.® 
Ignoring nonlinear effects until the PIT-crossing, par¬ 
ent growth is described by Qa = ^aQa, which means 
that PIT-crossing occurs at 


tpiT = 


I , Qpit 

— in ——^ 


(4.18) 


where Qq( 0 ) is the parent’s initial amplitude. 


F. Equilibrium solution 


Once the parent crosses the PIT and the daughters 
start growing, the three modes will continue interacting 
by exchanging energy. There can be two general out¬ 
comes from this process: (i) the system admits a stable 
equilibrium solution and all three modes reach satura¬ 
tion, or (ii) the parent’s growth cannot be halted by the 
daughters and all three modes grow, continuing to ex¬ 
change energy. 

Equations (4.11) admit an easy-to-obtain equilibrium 
solution. Expressing the complex amplitudes Q in terms 
of real amplitude and phase variables, we can introduce 
the variable transformation [52] 

Qa = (4.19a) 

ri 

(4.19b) 

Q 7 = (4-19c) 

Then, Eqs. (4.11) are written as 

= 7 a£a + sin:/j, (4.20a) 

£/3 = - £-i£a sin(4.20b) 

Er^ = — ea.£^ sin:/?, (4.20c) 


and 


(yj = cot (p 


^OL 

^OL 


ei 



+ Aw, 


(4.20d) 


where p = 'da—— 'd-y + Awt and 7 = 7 a + 7/3 + 77 - 
Setting the time derivatives to zero, we find the steady- 
state solution 


£« = 7/377 


e/3 = -777a 


= -7a 7/3 


1 + 1 A: 

7 


1 -b 


Aw 

7 

1 + 1^ 

7 


and 


cot p = 


Aw 

1 

7 


(4.21a) 

(4.21b) 

(4.21c) 

(4.21d) 


or, in terms of the original complex amplitudes, 

2 ‘ 

(4.22a) 
(4.22b) 


IQ |2 _ 7/377^/3^7 




1 + 1 ^ 

7 


\Qp? = - 

m^ = - 


T 7 T 0 ^7 ^0 

1 -h 

f Aw 



V 7 

lalpbabp 

1 + 

/^Aw 


V 7 


(4.22c) 


Note that, for [ 7,9 -|- 7 .^] ^ 7 a, the equilibrium ampli¬ 
tude (4.22a) of the unstable mode coincides with the PIT 
(4.17). 


G. Saturation conditions 

Such three-mode coupled systems, exhibiting a para¬ 
metric resonance instability, have been studied in the 
past [53, 54] for their significance in various fields, e.g. 
plasma physics [55, 56]. These studies show that certain 
conditions have to be met, in order for the system to 
approach saturation. 

Performing a linear stability analysis of Eqs. (4.20) 
(which is presented in Appendix C2), we find that the 
equilibrium solution (4.22) is stable if [52] 

\lp+li\>la (4.23) 


and 


® Note the importance of the mode frequency signs here: ifa;^aj-y < 
0, then and no parametric instability can occur. This is 

a result of the assumed resonance (4.7) between the parent and 


the daughters. If we perform the same analysis, for example, for 
mode /3 being the parent, then ps Uq, — cj-y, in which case 
< 0 is a necessary condition for parametric instability. 





























10 



FIG. 4. (a) A versus ((= (^ = ^j). The saturation condition (4.26) is satisfied inside the shaded area. The two asymptotes at 
( ~ 1.37 and A — 2( — 1 are also shown (dashed lines). A global minimum occurs at (1.77,3.73). (b) A versus (/s versus (j. 
The saturation condition (4.24) is satisfied inside the region that lies above the plotted surface. The thick line corresponds to 
the case where (p = 


3 I (C/3 + C 7 ~ 1) iC/3 — C 7 )^ + 2 (C/3 + C 7 ) + ^ 

+ {(C/3 + C 7 ~ 1) iCp ~ + (C/3 + C7)^ + 2 

- (C /3 + C 7 -1)' - 2 C/ 3 C 7 > 0 , 



(4.24) 


where C,a ,7 = — 7 / 3 , 7 / 70 ; which are the relative damp¬ 
ing rates of the daughters. To simplify the expression 
above, we set (^ = = (.y. Then, keeping in mind that 

Eq. (4.23) should also be true, it is reduced to 

C > « 1-37 (4.25) 

and 

^^>|J^|^(l-2Cr, (4.26) 

where A = Aw/ya. 

First, we notice that Eq. (4.25) imposes a stronger 
constraint on ( than Eq. (4.23). Second, we see from 
Eq. (4.26) that there is a lower limit on the detuning, 
which depends on (). This is illustrated in Fig. 4. 

If Eq. (4.23) is not satisfied, all solutions are un¬ 
bounded and the triplet’s amplitudes grow to infinity; 
the damping rate of the daughters needs to be larger than 
the driving rate of the parent, in order to stop its growth. 
The additional condition (4.24) [or, for 7 ^ = 77 , (4.25) 


and (4.26)] is more unintuitive; as shown by Refs. [53, 54], 
a number of interesting behaviors occur when it is not 
fulfilled, including limit cycles and chaotic motion. The 
amplitude evolution of a triplet satisfying the saturation 
conditions can be seen in Fig. 5. 

V. DISCUSSION 

The anticipated advent of gravitational-wave astron¬ 
omy will hopefully shed some light on the neutron star 
equation of state problem: should gravitational radiation 
from individual sources be observable, much information 
about the neutron star interior could be obtained. How¬ 
ever, gravitational-wave asteroseismology would have to 
deal with very weak signals, generated by stellar oscilla¬ 
tions. The fact that some of these oscillations are unsta¬ 
ble to the emission of gravitational radiation, due to the 
CFS mechanism presented in Sec. Ill, works to our ad¬ 
vantage: the amplitude of the mode will grow until such 
a point when nonlinear effects saturate the instability. 

Studies on the r-mode instability have shown that the 
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FIG. 5. Amplitude evolution of a coupled triplet that satis¬ 
fies the saturation conditions. Horizontal solid lines represent 
the saturation amplitudes of each mode. The dashed horizon¬ 
tal line shows the position of the triplet’s PIT, whereas the 
dashed vertical line denotes the PIT-crossing time. At that 
point the parent (mode a) crosses the PIT and the daughters 
(modes /3 and 7 ), which were damped until that point, start 
to grow. Then, the amplitudes oscillate and finally converge 
(albeit very slowly in this example) around their equilibrium 
values (the parent’s equilibrium coincides with the PIT in 
this example). In this graph, we show the triplet with the 
lowest PIT, in a polytropic model with P = 2 and Pi = 2.1, 
for n = Hk and T = 5 x 10® K. Mode a is the |/-mode, 
mode P is the ~ 4 /-mode, and mode 7 is the ygs-mode (where 
the notation has been used). The growth/damping rates 
are 7 ^ = 2.7 x 10“®rads“^, 7/3 = —1.0rads“^, and 7 -, = 
— 1.4 X 10~®rads“^, and the detuning is Ao; = 14.1rads“^. 
The value of |Q| depends on the mode normalization choice 
as \Q\ — -JEjnode/Eunit (cf. Sec. IV C); here, we chose 

Sunit = Mc^. 


and two selection rules for their orders m and degrees 1. 
Although any triplet can be part of this network, we are 
obviously interested in the case where one of the partici¬ 
pating inodes is the unstable /-mode. Then, the coupled 
triplet is said to be parametrically resonant and can lead 
to a parametric instability, if the unstable (parent) mode 
crosses the so-called parametric instability threshold. At 
that point, the other two (daughter) modes start grow¬ 
ing. The system reaches saturation if certain conditions 
are satisfied for the modes’ growth/damping rates, and 
their frequency mismatch. 

In this paper, we have focused on polar modes, like /-, 
p-, and g-modes. However, all the formulas presented in 
Sec. IV are also applicable to axial modes. It is only in 
Appendix B where we assume that all three modes are 
polar, and find an expression for the zeroth-order compo¬ 
nent of the coupling coefficient. Results from the appli¬ 
cation of the formulation above to Newtonian, polytropic 
stars will be presented in a subsequent paper. 
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APPENDIX A: DERIVATION OF THE 
EQUATIONS OF MOTION 

1. The quadratic equation of motion 


saturation levels will make detection very difficult. In the 
most optimistic cases, the signal may be detectable with 
Advanced LIGO from within the local galaxy group [23]. 
As far as the /-mode instability is concerned, reasonably 
high saturation levels make the signal from a nascent star 
definitely detectable with the Einstein Telescope (in some 
cases even with Advanced LIGO) for sources in the Virgo 
cluster [32]. 

Estimating the saturation amplitudes for the r- and 
/-mode instabilities is also important for another reason: 
their values affect the evolution of the star inside the in¬ 
stability area. A newborn star, for which both instabili¬ 
ties can be significant, will enter the instability window, 
which it will traverse at approximately constant angu¬ 
lar velocity, until it reaches thermal equilibrium; then, at 
approximately constant temperature, the star will spin 
down due to the emission of gravitational radiation, as 
well as magnetic braking, until it exits the window. The 
saturation amplitude affects the duration of these phases, 
thus the time which the star spends inside the instability 
area. 

By taking quadratic perturbations into account, cou¬ 
pled three-mode networks are formed throughout the 
star. These triplets have to satisfy an internal resonance 


The derivation of the quadratic equation of motion 
(4.3) can be performed in the same way as the derivation 
of the linear equation of motion (2.10), except that now 
we also want to retain second-order perturbative terms. 

Following Ref. [52], we will use the velocity v, instead 
of the Lagrangian displacement to describe the pertur¬ 
bation. As mentioned in Sec. II, the background velocity 
is zero (because we are working in the corotating frame), 
so V = 6v = ^. Differentiating Eq. (2.3) with respect 
to time and imposing perturbations on the equilibrium 
state, we obtain the equation of motion for the velocity, 
namely 

V + B{v)+C{v) + = 0, (Al) 


where 


B{v) = 2Q X V (A2) 


and 





VpdSip / 9(51$ \ 

dt \ dt )' 


(A3) 


with i5i denoting first-order and ^2 second-order Eulerian 
perturbations. Afv represents the quadratic terms, which 
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are explicitly written as 


and 




(v ■ 'V)v + 


where 




dip 


C) 


V 62 P 

p 


+ (5i - V(5ip 


Vp + V(52$ 


(A4) 


1^1 - =-^ and (52 - =-^ 


d 2 P , {5ipY 


d 62 P 

where 


-( 7 ; • V) V + [(I • v)(pri) +prixv • |] v • i;, 

(All) 


X = Ti + 


/ainTi 
\ 91np 




2. The amplitude equation of motion 


It should be noted that AT, which appears in Eq. (4.3), 
is related to simply by = dAf/dt. 

Perturbing the continuity equation (2.2), we get 

^ = -py-v-{v V)p (A5) 

and 

^ = -SipV -v-iv V) 6 ip, (A6) 

for first- and second-order terms, respectively. Accord¬ 
ingly, the perturbed Poisson equation (2.4) gives 

= AnGdip and ^^ 62 ^ = 47rG'(52p, 


In order to obtain the equation of motion for the am¬ 
plitude (4.4), we have to replace v in Eq. (Al) with the 

expansion (4.2). Note that this expansion implies that 

E = 0 (A12) 

OL 

and 

O' 

= 0. (A13) 


whose (time-differentiated) solutions are 


and 



|r - r'\ 


(A7) 


Making use of the eigenvalue equation (2.13), the orthog¬ 
onality condition (2.20), as well as Eqs. (A12) and (A13), 
we get 


Qo 


■iuJaQa = J-{^a,Af^)e 
On 


(A14) 


9^2 4* 
dt 



Vr' ■ {dipv) ^3^> 


(AS) 


Finally, perturbation of the equation of state p = 
p{p, p) to second order gives 



(A9) 


where Pi is defined by Eq. (2.9). Here we have assumed 
that Ap = 0, i.e. the composition is frozen (if p corre¬ 
sponds to the composition) and/or the star is isentropic 
(if p denotes entropy). Also, we have used Lagrangian 
perturbations, which, to second order, are related to Eu- 
lerian by 

A/ = 5if + {^-V)f + 62 f + {^-V)5if + \^-[^ ■ V (V/)]. 


Using this, we obtain from Eq. (A9) 


ddip 


— {v ■ V)p — pPiV ■ V 


(AlO) 


It is easily seen that Eq. (A14) is obtained by differenti¬ 
ating Eq. (4.4) with respect to time. By further replacing 
the expansion (4.2) in Afv one gets 

Qn+ioJnQa = 

“ /3.7 

(AI5) 

where 

Fap-y = ^—{$,n,Af.v{^P,$,^)) (AI6) 

is the coupling coefficient (a bar over an index means that 
the corresponding mode eigenfunction in Af-u has to be 
complex conjugated and its frequency sign reversed). 

As mentioned in Sec. IV B, not all terms in Eq. (AI5) 
play an equally important role in the amplitude evolu¬ 
tion. As shown in Appendix C 1, a resonance condition 
between the modes is necessary for the dynamics of the 
system to be significantly affected by quadratic terms. 
Assuming a resonance of the form Wq, = uip + uj^ + Aw, 


dt 
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where Aw is a small detuning, one can omit rapidly vary¬ 
ing terms in Eq. (A15). Then, choosing a mode triplet 
which satisfies the resonance condition, we get 


Qa j (A17a) 

Oq^ 

Qp + itupQp = (A17b) 

0/3 

(A17c) 

If such a resonance exists, it can be shown that 
iuJaFap-y = i{uja — Auj)Fap-y, where Fap-y is given by 
Eq. (4.6). So, ignoring the detuning, Fap^ « Fap~f, which 
also implies that Q is negligible, because only then we can 
retrieve the equivalent system (4.8). 

Setting H = F^p^ = Fp^a = F^^p (cf. Appendix B) 
and introducing growth/damping rates for the modes. 


Eqs. (A17) become 


jqj 

Qc = 7 aQa + -f-QpQjC-^^^^ 
Oq/ 

, (A18a) 

jqj 

Qp = ipQp + 

(A18b) 

jqj 

Qj = 77O7 + 

(A18c) 

which coincide with Eqs. (4.11). 


APPENDIX B: THE COUPLING COEFFICIENT 

Proceeding with the evaluation of Eq. 
equations from Appendix A 1 , we find an 
for the coupling coefficient, which is [52] 

(A16), using 
explicit form 


FaP'Y i^^p^aP'y 4“ , (^^1) 

COa 

where 




J j (^/3 • ^ 7 ) + 1/3 X (V X X {\7 X ^p) ] 


- -[v • (pC/ 3 ) V • Vp + pEiV • 1 ^,) + V • {p^j)V i^p -Vp + pTiV ■ |; 3 )] 


+ V-(p£/3)V-(p|^)^- 


1/3 -V 


V• {p^,) 


Vp - GpV 


Vr' • [4/3 V • (^^ 7 )] 


r — r' 


+ V [1/3 • V (1^ • Vp + pTiV • + (V • ip) ■ V (pTi) + pEix (V • ip) (V • l^) ] • CP?r. (B2) 


The expressions for Fp^a and F^^p are obtained from 
Eq. (Bl), keeping in mind that a bar over an index 
means that the corresponding mode eigenfunction has to 
be complex conjugated and the corresponding frequency 
has to change sign. 

As pointed out by Ref. [18], the expression above for 
the coupling coefficient is identical for both nonrotating 
and rotating stars. This of course does not make the ac¬ 
tual value of the coupling coefficient the same for both 
cases. If rotation is included, the eigenfrequencies, the 
eigenfunctions, and the equilibrium quantities are all af¬ 
fected (cf. Sec. IIB). 


ties [ 1 ] 


X = r/R, 


Vi = 

2/3 = 


r 

54) 


gr 


Cl = 


Q 


i 1 


U}=Ujly^GM/W, 
. 2 C /1 


2/2 = CiW 

2/4 = 


r 

1 

9 dr ’ 



U = 

dlnr 

Mr ’ 

1 dlnp 

A* = 

1 dlnp 

dlnp 

Ti d In r ’ 

Ti dlnr 

dlnr 


We now assume that ^ takes the form (2.14), namely, 
it describes the eigenfunction of a polar mode in the non¬ 
rotating limit. We also define the dimensionless quanti- 


where g = GMr /is the local gravitational acceleration 
and Mr = 47rpr^dr. Then, after cumbersome calcula¬ 
tions, the coupling coefficient takes the form [52] 
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n = 


n 


GM/R^ 


^z. 


a/S-y 


f . / QCk 

/ ^ {^*yi,k + VgZk) ( wk'Wk"yi,k'yi,k" + -^^3 

•^0 I fe ^ 


■ 2 / 2 ,fc' 2 / 2 ,fe" 


Cl 

Cl 

>1 


dlnFi 




Vg + U - 4:- Ci'^ml\Y\_ 2 / 1 ,fc - X! ^kyi,k'yi,k" + ^ 2 / 4 .fe 2 /i,fe' 2 /i,fc" 

V. k / k k k 


H—2" X! iGCkyi,k'yi,k" + QGk'yi,k'Zk" + QCk"yi,k"Zk') \ pR^x'^Ax. 


(B3) 


In the expression above, the index k successively takes 
one of the values {a, /3, 7 ), whereas the indices k' and k” 
take the values that come next and after next, respec¬ 
tively (for example, for k = a, k' = 13 and k” = 7 ). The 
rest of the quantities are defined as 


Zk = 2/2,fc — 2/3,fc5 


'^k 


uJk r k = a 

~ 10 ^ 1 a 

-UJk k = p,^, 


QGk 


—Afc -|- Ak' + Ak" 

‘IWk' ZUk^’ 


GGk 


AkZUk {Ak' Ak") {zUk' ZUk") 

‘2zU kZXJ k'ZXJ k" 


with Ak = Ik (Ik + !)■ Also, 


A 


g 


dlnTi / SlnTiN 

dlnr ® ainp ) 


Finally, 


Z 


Oi^y — 


Y*YpYy sin OAOdcf), 


where Yk = . 

Equation (B3) is invariant to the transformations 

Ya < Zp-, yi^a ^ yi,Pi Zy / Yy , UJj / UJj 

and 

Aq: Yy, yi,a ^ 2/1,7; Ay 3 ^ Yg , UJp ^ ^/ 3 ; 

which proves that Fap^ = Fp^a = F-yap = F-- 

The expression above is the zeroth-order component of 
the coupling coefficient, namely, all quantities are evalu¬ 
ated in the nonrotating limit. A more general expression 
could be found if we had replaced the rotationally cor¬ 
rected eigenfunctions in Eq. (Bl), but this would signifi¬ 
cantly complicate the calculation. 


R has units of energy; the normalization in Eq. (B3) 
is useful when all quantities in the amplitude equations 
of motion (4.11) [or (A18)] are normalized accordingly. 
Defining a dimensionless time r = ty/GM/R^ and a di¬ 
mensionless frequency w = ujjy/cM/lF, the equations 
of motion are written 


Q'c = aMc + (B4a) 

^OL 

jqj 

Q'p = ipQp + (B 4 b) 

bp 

iR 

g; = 77 Q 7 + yQcQpF^‘^\ (B4c) 

where 7 = "f/^^GM/R^, b — b/y/GM/R^ and the prime 

denotes differentiation with respect to t. 


APPENDIX C: STUDY OF A THREE-MODE 
NETWORK WITH QUADRATIC 
NONLINEARITIES 


1. The multiscale method 


Let us assume that we have an ordinary differential 
equation which includes a small parameter e. We write 
the solution to this equation in the form of an asymptotic 
series, in the sense that 


00 

y{t) '^yn{t)e^- 

n—0 

In the beginning of the evolution, when t is small, low- 
order terms dominate the solution. However, as t grows 
bigger, the contribution of higher-order terms cannot be 
neglected. These terms are usually called secular terms, 
because their effects become important (compared to 
low-order terms) at later stages of the evolution. This 
behavior appears, for example, in a damped harmonic 
oscillator, where the zeroth-order solution is simply an 


















15 


undamped harmonic oscillation, with the damping effects 
occurring at higher orders. 

The multiscale method (cf. for instance, Ref. [48]) is 
a way to capture such higher-order effects from secular 
terms and make them appear in the low-order terms. 
As a result, the low-order approximation of the solution 
would be valid on secular time scales. 

We define the time scales Tn = e”t and rewrite the 
asymptotic solution, so that 

OO 

n=0 


Replacing the solutions in Eqs. (Cl) and distinguishing 
between 0 {e) and C>(e^) terms, we get 


dTo 




= 0 , 
= 0 , 
= 0 , 


In other words, we let the terms of the series depend on 
more than one time scale. As we will see, this allows 
us to “eliminate” secular effects from higher-order terms, 
thus preventing these terms from becoming signihcant. 

We are going to use this method, in order to study 
Eqs. (4.11). First, we remove the exponential time de¬ 
pendence by setting Ck = Qk exp{iuikt) (k = a, /3, 7 ) and 
the equations of motion are written as 

iV 

Ca — iiOaCa = ^aCa + (Cla) 

iV 

Cp - iujpC0 = -fpCp + (Cib) 

i'h 

+ —CaC}. (Clc) 

Now, we seek solutions of the form 


and 


dC[ 


( 1 ) 


da 


( 2 ) 




dT. 






( 1 ) 


-.( 2 ) 


dTy 




dc: 


( 1 ) 


dTi 


0 

( 2 ) 
7 

dTn 


da 


respectively, where we also set 7 ^ = £%, so that damping 
and nonlinear terms appear in the same order. 


The first-order equations have simple solutions of the 
form 


Ck = eC'W(To, Ti) + e^cj^\To, T,) + ..., 
where Tq = t and Ti = et. Time derivatives then become 


cfHTo,ri)=Afc(ri)e“'=^“, (C2) 


JL 


which we substitute to the second-order equations, to get 


dC, 


( 2 ) 


dTo 

dCf 


da 


( 2 ) 


dTo 


- iuJaCi^^ 

- iupCf^ 


i^a-^OL 

( 77^7 


dAa\ 

drj 

cL4^\ 

drj 

dA.^\ 


Uot 

j'-U 

bp 

jnj 


(C3a) 

(C3b) 

(C3c) 


As we mentioned earlier, the whole point of the multi¬ 
scale method is to transfer long-term effects from higher- 
order terms to low-order terms. In this case, we want 
to prevent the second-order terms of the solution, 7 
from growing and becoming important. To accomplish 
this, we have to eliminate the so-called secular terms. 
In the case of Eqs. (C3), terms that include the factor 


exp {iojkTo) have to vanish, because they produce secular 
terms, causing the solution to grow in time. 
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a. The nonresonant case 


If there is no resonance of the form « W /3 + ui-y be¬ 
tween the modes, then the conditions for the elimination 
of secular terms from Eqs. (C3) are 


dAfc 

dTi 


IkA-k: 


or 


Ak = ake"‘'^^, 

which makes the first-order solutions (C2) 

Ck = eCf) + 0 {e^) = + 0 {e^), 

or, in terms of the original variables Qk, 

Qfc =eafee^''*+0(e2). (C4) 

Equation (C4) shows that, if there is no resonance be¬ 
tween the modes, their amplitudes grow or decrease with 
time, depending on the sign of jk ■ 


b. The resonant case 


If a resonance of the form uJa = + oj-y + Aw exists 

(Aw being a small detuning), then the second terms on 
the right-hand sides of Eqs. (C3) also contribute in the 
production of secular terms in the solution. Then, the 
secular-term elimination conditions become 


dAa 

dTi 


%A, 


+ 7 — 

Ury 


— iAc 2 jTi 


dA^ 

dTi 

dA^ 

m 


= ipA^ + 

= %A^ + '^Ao,Ay^^^\ 


(C5a) 

(C5b) 

(C5c) 


where we set Aw = eAw. From Eqs. (C5) we obtain 
our original system (4.II), whose study is presented in 
Secs. IVE-IVG. 


2. Linear stability analysis 


Having used the variable transformation (4.19) to the 
equations of motion (4.11), we obtain Eqs. (4.20), namely 


= la£a +epe~fSinip, 
£/3 = 7/3£/3 - e7£asin(/j, 
£7 = — £q ,£/3 sin(/J, 


and 


If = cot (f 


£oi 

£q 


ei 

£/3 



-I- Aw, 


where ip = 'da — 'dp — '&-y + Aujt and 7 = 7 a + 7/3 + 77 - 
We linearize these equations by imposing small pertur¬ 
bations around their equilibrium solutions (4.21). Denot¬ 
ing these perturbations by <5 (not to be confused with a 
Eulerian perturbation), we get [52] 


and 


d / Sec 
dt \ £a 

d / Sep 
dt V ep 

d / ^£7 

7 


= -7a 


= -7/3 


= -77 


fee 

Sn/ 


Sep 

ep 




SSf^ Ssp 




6 e^ 


kS(p 


■ kS(p 
tvSip 


dST Y- r , X 

= k} Tk - \-l 0 ip, 

k 


dt 


(C 6 a) 

(C 6 b) 

(C 6 c) 

(C 6 d) 


where k = Aw /7 Fk = 2jk — 7 , with the index k 
successively taking the values (a,/ 1 , 7 ). 

The matrix of the linear system (C 6 ) is 


/ 7a -7a -7a \ 

-7/3 7/3 -7/3 -K7/3 

-77 -77 77 -«77 I ’ 

V nFa nFp kF^ 7 / 

with the help of which we can find the system’s character¬ 
istic polynomial, via the relation \A — A/| = 0, where A 
are the eigenvalues of A and I is the identity matrix. The 
polynomial has the form A^ -I- oi A^ -I- a 2 A^ -I- a^X -I- 04 = 0, 
where 

ai =- 27 , 02 = 7^ (1 + ^^) - Iklk ', 

k 

03 = 4 (1 + 3 /^^) 04 = -4 (1 -h 7 

k k 


with the index k' taking the value that comes after Fs 
value (e.g., if fc = a, k' = (3). 

Now, we can use the Routh-Hurwitz stability criteria 
(cf. for instance. Ref. [57]), in order to determine the 
behavior of the system. First, we construct the Routh- 
Hurwitz matrix, using the polynomial coefheients, which 
is 


M = 


/ oi 1 0 0 

03 02 Oi 1 
0 04 03 02 

Vo 0 0 04 


Then, the stability criteria are given by 

lEi = oi > 0 , 

Oi 1 

03 02 

oi 1 0 

03 O2 Oi 
0 04 03 


IE 2 = 


IE 3 = 


— 0402 — 03 > 0, 

= 0311^2 — 0^04 > 0, 


(C7) 

(C 8 ) 

(C9) 
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and 

W 4 = \M\ = a^Ws > 0. (CIO) 

Since < 0, it can be easily shown that the second 
and fourth criteria are redundant and follow from the 
other ones. Indeed, if Wi > 0 then 04 is also positive, 


which, combined with W 3 > 0 , makes the fourth criterion 
true. Also, W 3 > 0 yields W 2 > 0 ^ 04 / 03 , but since 
03 > 0 , the second criterion is also true. 

So, finally, from the first and third criteria we obtain 
the stability conditions (4.23) and (4.24), which are fur¬ 
ther studied in Sec. IV G. 
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